#!/usr/bin/env python
# 下届
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import ConnectionPatch



import matplotlib.pyplot as plt
from scipy import integrate, inf
from matplotlib.patches import ConnectionPatch
x = np.linspace(0.01, 5, 15)


def to_be_integrate(t, _x):     # 被积函数
    return t * np.exp(-(t**2)) * np.exp(-2*t*_x)


x2 = x * x
e2 = np.exp(-x2)
sqrtPi = np.sqrt(np.pi)
B1 = [(1 - 2 * integrate.quad(to_be_integrate, 0, inf, args=(i, ))[0]) for i in x]
B2 = e2/(np.sqrt(np.pi)*x)
erfc = B1 * B2

C1 = -(1+2*x)
# r5 = B2 * (1 - ((2 - (3+2*x)*np.exp(C1)) / ((1+2*x)**2)))
# r6 = B2
# r7 = e2
# r13 = (2 * (x2+2) * np.exp(-(x2/2))) / (np.sqrt(np.pi) * x * (x2+3))
# r27 = 2*np.sqrt(2) * (e2/50 + np.exp(-(x2/2))/(2*(x+1)))


# 下界
r24 = e2 * (1 / (sqrtPi*x)) * (1 - ((-C1) + (3+2*x)*np.exp(C1))/((1+2*x)**2))
r26 = 2 * ((1/12)*e2 + (np.exp(-(x*x/2)) / (np.sqrt(2*np.pi)*(x+1))))
r27 = 2 * e2 / (sqrtPi * (x + np.sqrt(x*x+2)))
fx = np.sqrt(x**4 + 6*(x*x) + 1)
r29_1 = np.exp((fx+x*x-1)/(-4))
r29_2 = np.sqrt((np.e*(fx+x*x+1))/np.pi)
r29_3 = fx + x*x + 3
r29 = 2*r29_1*r29_2/r29_3

plt.figure(figsize=(16, 8), dpi=98)
p1 = plt.subplot(121, aspect=5 / 2.5)
p2 = plt.subplot(122, aspect=0.5 / 0.05)

line_width = 2
p1.semilogy(x, erfc, '-0', lw=line_width, label='erfc')
p1.semilogy(x, r24, '-1', lw=line_width, label='we proposed lower bound')
p1.semilogy(x, r27, '-3', lw=line_width, label='the complementary error function')
p1.semilogy(x, r26, '-4', lw=line_width, label='very simple tight bounds')
p1.semilogy(x, r29, '-v', lw=line_width, label='A Tight Lower Bound on the Guassian Q-Function')

p2.semilogy(x, erfc, '-0', lw=line_width, label='erfc')
p2.semilogy(x, r24, '-1', lw=line_width, label='we proposed lower bound')
p2.semilogy(x, r27, '-3', lw=line_width, label='the complementary error function')
p2.semilogy(x, r26, '-4', lw=line_width, label='very simple tight bounds')
p2.semilogy(x, r29, '-v', lw=line_width, label='A Tight Lower Bound on the Guassian Q-Function')

p1.axis([0.0, 5.01, 1 / (10**10), 1])
p1.grid(True)
p1.legend()

tx0 = 3.2
tx1 = 3.25
ty0 = 8 / (10**6)
ty1 = 4 / (10**6)

p2.axis([tx0, tx1, ty1, ty0])
p2.grid(True)
p2.legend()

sx = [tx0, tx1, tx1, tx0, tx0]
sy = [ty0, ty0, ty1, ty1, ty0]
p1.plot(sx, sy, "purple")

xy = (tx1, ty1)
xy2 = (tx0, ty1)
con = ConnectionPatch(xyA=xy2, xyB=xy, coordsA="data", coordsB="data",
                      axesA=p2, axesB=p1)
p2.add_artist(con)

xy = (tx1, ty0)
xy2 = (tx0, ty0)
con = ConnectionPatch(xyA=xy2, xyB=xy, coordsA="data", coordsB="data",
                      axesA=p2, axesB=p1)
p2.add_artist(con)


# plt.ylim(np.log10(1/(10**10)), 0)
# plt.legend()

# plt.savefig('compare.eps', format='eps')
# plt.savefig('compare.eps', dpi=60, format='eps')

plt.show()
# out_fig = plt.gcf()
# out_fig.savefig('compare.eps', dpi=600, format='eps')

